#identify modules

# Rscript step2.R --exprFile "../test_data/datExpr.csv" --network "../1_construct_network/network.RData" --output "../2_identify_modules"

library(optparse)


option_list = list(
  make_option("--exprFile", type="character", default=NULL,
              help="expression file path"),
  make_option("--network", type="character", default=NULL,
              help="network file path"),
  make_option("--selectGeneNum", type="integer", default=10,
              help="select gene number [default= %default]"),
  make_option("--output", type="character", default="out.txt",
              help="output path [default= %default]")
)

# 解析命令行参数
opt_parser = OptionParser(option_list=option_list, add_help_option=TRUE)
opts = parse_args(opt_parser)

datExpr = read.csv(opts$exprFile, stringsAsFactors = FALSE, row.names = 1);

library(WGCNA)

load(opts$network)
nGenes = length(net$goodGenes)

moduleColors = labels2colors(net$colors)

# Visualizatioin network modules
png(
  filename = paste(opts$output, "network_module.png", sep="/"),
  width = 13,
  height = 9,
  units = "in",
  res = 300
)
plotDendroAndColors(net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)

dev.off()


# Visualizatioin network heatmap
png(
  filename = paste(opts$output, "network_heatmap.png", sep="/"),
  width = 13,
  height = 9,
  units = "in",
  res = 300
)
dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = 6);

# 随机选择部分基因
nSelect = ifelse(opts$selectGeneNum > nGenes, ceiling(nGenes / 10), opts$selectGeneNum)
set.seed(10);
select = sample(nGenes, size = nSelect);
selectTOM = dissTOM[select, select];
selectTree = hclust(as.dist(selectTOM), method = "average")
selectColors = moduleColors[select];
plotDiss = selectTOM^7;
diag(plotDiss) = NA;
TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")

dev.off()
